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Motivated by questions about the QCD deconfining phase transition, we studied in two previous 
papers Model A (Glauber) dynamics of 2D and 3D Potts models, focusing on structure factor 
evolution under heating (heating in the gauge theory notation, i.e., cooling of the spin systems). 
In the present paper we set for 3D Potts models (Ising and 3-state) the scale of the dynamical 
effects by comparing to equilibrium results at first and second order phase transition temperatures, 
obtained by re-weighting from a multicanonical ensemble. Our finding is that the dynamics entirely 
overwhelms the critical and non-critical equilibrium effects. 

In the second half of the paper we extend our results by investigating the Glauber dynamics of 
pure SU(3) lattice gauge on N T N„ lattices directly under heating quenches from the confined into 
the deconfined regime. The exponential growth factors of the initial response are calculated, which 
give Debye screening mass estimates. The quench leads to competing vacuum domains of distinct Z3 
triality, which delay equilibration of pure gauge theory forever, while their role in full QCD remains 
a subtle question. As in spin systems we find for pure SU(3) gauge theory a dynamical growth of 
structure factors, reaching maxima which scale approximately with the volume of the system, before 
settling down to equilibrium. Their influence on various observables is studied and different lattice 
sizes are simulated to illustrate an approach to a finite volume continuum limit. Strong correlations 
are found during the dynamical process, but not in the deconfined phase at equilibrium. 

PACS numbers: PACS: 05.50.+q, 11.15.Ha, 12.38.Gc, 25.75.-q, 25.75.Nq 



I. INTRODUCTION 

In investigations of the QCD deconfining phase tran- 
sition (or crossover) by means of heavy ion experiments, 
one ought to be concerned about non-equilibrium effects 
due to the rapid heating of the system. With this in 
mind we have investigated in previous papers Q, Q] the 
Model A (Glauber) dynamics of 2D and 3D Potts 
models. Model A dynamics includes all diffusive stochas- 
tic local updating schemes (Metropolis, hcatbath, etc.) 
and not only the process introduced in In 3D Potts 
models spins provide degrees of freedom, which mimic 
Polyakov loops effectively p|, while in 2D analytical re- 
sults [f| allow to check on the accuracy of the employed 
numerical methods. For other approaches to simulate 
non-equilibrium quantum fields see Ref. 0. 

The QCD high temperature vacuum is characterized 
by ordered Polyakov loops, which are similar to spins in 
the low temperature phase of the 3D 3-state Potts model. 
We model heating by a quench from the disordered into 
the ordered phase, which thus corresponds to a cooling 
quench in the analogue spin model. Time evolution after 
the quench leads to vacuum domains of distinct triality 
under the Z 3 center of the SU(3) gauge group. It ap- 
pears that these competing domains are the underlying 
cause for the explosive growth of structure factors Fi(t), 
which we encounter in the time evolution after a heating 
quench. We use the term spinodal decomposition loosely 
to denote generically such a time period of globally un- 
stable behavior. 

Relaxation of the system at its new temperature be- 
comes only feasible after each structure factor has over- 



come its maximum value. While the maximum value of 
the structure factor diverges with lattice size, its initial 
and final equilibrium values are finite in the normaliza- 
tion chosen in the paper. The time (measured in updates 
per degree of freedom) for reaching the maximum di- 
verges with lattice size unless the underlying order-order 
symmetry is broken. Once the system has equilibrated at 
high temperature, the subsequent temperature fall-off is 
driven by spatial lattice expansion and the system stays 
in quasi-equilibrium during this period. So one has dif- 
ferent time scales under heating and cooling 

The early time evolution of SU(3) gauge theory af- 
ter the quench is well described by stochastic equations, 
which follow from dynamical generalizations of equilib- 
rium Landau- Ginzburg effective action models. We cal- 
culate the exponential growth factor of this linear ap- 
proximation and use a phenomenological model 9] to es- 
timate the Debye screening mass for two temperatures 
above the deconfining T c . 

Finally we compare measurements of Polyakov loop 
correlations, gluonic energy densities and pressures 
around structure function maxima with their equili- 
brated values in the deconfined region at high tempera- 
tures. These measurements are of interest for a scenario 
in which the heating process turns back to cooling before 
actually reaching the equilibrium side of the structure 
factor maxima. In the conclusions we continue this dis- 
cussion. 

In the next section we introduce our notations and 
some preliminaries. Section IlIII deals with Potts models. 
First equilibrium properties of structure functions are es- 
tablished by means of multicanonical simulations of the 
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3D Ising and 3-state Potts model. Subsequently their dy- 
namical evolution after a quench is investigated, extend- 
ing previous results. In section HVI we present our simu- 
lations of pure SU(3) lattice gauge theory. Some SU(3) 
data were already reported at the 2004 APS DPF confer- 
ence [lOj. As these simulations are very CPU time con- 
suming it took over one more year to collect the present 
statistics. Summary and conclusions are given in the final 
section Ivl 



II. NOTATION AND PRELIMINARIES 

We summarize our basic notations and concepts in this 
section. 

A. Models 

We simulate g-state Potts models with the energy func- 
tion 



B. Structure Factors 



E = 2 E(V*» 



(1) 



where the sum is over nearest neighbors of a hypercubic 
lattice in D dimensions. The spins qi of the system take 
on the values <?, = 0, . . . , q — 1. The factor of two and 
the term 1/q is introduced to match for q = 2 with Ising 
model conventions Simulations are carried out with 
the Boltzmann factor exp(— (3E). 

The Wilson action for pure SU(3) non-Abelian Eu- 
clidean lattice gauge theory is 



where U n 



7l,fJ,lS 



(2) 



L/t denotes the prod- 



uct of the SU(3) link matrices in the fundamental rep- 
resentation around a plaquette and the sum runs over 
all plaquettes. Simulations are done with the Boltzmann 
factor exp(5' J 4). 

The Markov chain Monte Carlo (MC) process pro- 
vides model A (Glauber) dynamics in the classification of 
Ref. |a| . For Potts models we use the heatbath algorithm 
of [ill and for SU(3) gauge theory the Cabibbo-Marinari 
|l2fl heatbath algorithm and its improvements of Ref. 0] 
(no over-relaxation, to stay in the universality class of 
Glauber dynamics). In each case a time step is a sweep of 
systematic updating through the lattice, which touches 
each degree of freedom once. With small statistics we 
have checked that updating in random order gives sim- 
ilar results up to a slowing down of the evolution speed 
by a constant factor. This is expected as in equilibrium 
simulations random updating has larger autocorrelations 
than systematic updating |l l| . For our equilibrium sim- 
ulations of Potts models we used a multicanonical 0] 
Metropolis algorithm. 



Consider two-point correlation functions defined by 



where i denotes spatial coordinates. Periodic boundary 
conditions are used and the subscript L on the left-hand 
side reminds us that the average is taken over the spatial 
lattice. For gauge systems we deal with fluctuations of 
the Polyakov loop, for analogue spin systems with fluc- 
tuations of the magnetization. 

The finite volume continuum limit of @ is achieved 
by lattice spacing a — > and N a — > oo with the physical 
length of the box L = aN n = const. This means that 



a - ,V cr - 



transforms into 



(u(0)u t (i?)) L = -3 / d z ru{r)u\r + R), 



(4) 



(5) 



with f=ai, R = aj, u(r) = ito(z), and so on. We define 
the structure function F(p) as Fourier transformation of 
the two-point correlation function JSJ: 



F(p) 



[u(0)u\R)) L e l P R d 3 R. 



Periodic boundary conditions imply: 

_ k 2vr 
P = ~ = t"i 



(G) 



(7) 



where n is an integer vector (0, 0, 0), (0, 0, 1), and so on. 
The discretized version of © is 



F(p) = J> 3 (u (0)ul(j)) L e i% l 



(8) 



Using the definition and shifting the j summation one 
arrives (after straightforward algebra) at the expression 



E 



e- lkt u (i] 



(9) 



where we may rewrite the product in the exponent as 

2tt ^ 



api 



k 1 



(10) 



As we let the system evolve after a quench uo(i) becomes 
time-dependent: ua(i,t). The time t corresponds to the 
dynamical process, i.e., in our case the Markov chain 
model A dynamics. We consider an ensemble of systems 
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(replica) and dynamical observables are calculated as en- 
semble averages denoted by (...). The time-dependent 
structure functions averaged over replicas are: 



F ff (t) = (F(p,t)) 



(11) 



During our simulations they are averaged over rotation- 
ally equivalent momenta and the notation 



is used for the structure function at momentum 



2?r 



P = 



(12) 



(13) 



where \n\ — defines i. The Fi are called struc- 
ture function modes or structure factors (SFs). We 
recorded the following modes (including the permuta- 
tions) nr. (1,0,0), n 2 : (1,1,0), n 3 : (1,1,1), n 4 : (2,0,0), 
n 5 : (2,1,0), ne: (2,1,1), n 7 : (2,2,0), n 8 : (2,2,1) and 
(3,0,0), ng: (3,1,0). Note that there is an accidental 
degeneracy in length for rig. We measured also higher 
modes, in some cases up to n^. They exhibit the same 
behavior as the lower modes, but the data are far more 
noisy, so that we abstain from reporting these results. 
A difference to the normalization of Q, Q is that in the 
present paper we average over the permuted momenta 
instead of just summing them up. For instance, for the 
F\ SF the difference is a multiplicative factor of three. 



TABLE I: Statistics and SF maxima Fj max at f3 m from our 
equilibrium simulations of the 3D Ising model on L 3 lattices. 



L 


sweeps 






cycles 


20 


32 x 5 ■ 10 4 


17.00 (26) 


0.219874 


83 


30 


32 x 2 ■ 10 5 


37.97 (70) 


0.220825 


89 


44 


32 x 6 ■ 10 5 


78.4(1.6) 


0.221146 


70 


56 


32 x 1 ■ 10 6 


129.6(2.7) 


0.221345 


45 


66 


32 x 1.6 ■ 10 6 


175.6(4.9) 


0.221387 


43 


80 


64 x 2 ■ 10 6 x 3 


257.4 (2.4) 


0.221462 


65 + 72 + 67 



multicanonical parameters were extracted by finite size 
(FS) extrapolation from smaller to larger system, which 
is an efficient way when the FS behavior is controllable. 

The advantage of using multicanonical simulations is 
that accurate values of the SF peaks can be determined 
from one data set. Re- weighting of a canonical simulation 
1 1 ol ] allows accurate determination of the maxima of one 
quantity, but on finite lattices the maxima of different ob- 
servables are too far apart to be within the re-weighting 
range of one canonical simulation. We find it convenient 
to have the entire range of interest covered in one sim- 
ulation. In particular equilibration of the configurations 
around the transition and in the ordered phase is then 
secured due to frequent excursions into the discorded re- 
gion all the way to j3 = 0. 



III. POTTS MODELS 

For the analogue spin models the lattice spacing a can- 
not be varied. We set a = 1, so that the distinction be- 
tween L and 7V CT , p and k becomes superfluous. We use 
L and k in the following. The normalization of the SFs 
differs from our previous work QLQ. It is chosen so that 
they approach constant values in the infinite volume limit 
of equilibrium simulations of spin systems at non-critical 
temperatures. This follows from the fact that the ran- 
dom fluctuations in l|12|l are of order \f~V = \flfi. At a 
critical temperature of a second order phase transition a 
divergence of the SFs is then encountered as we illustrate 
for the 3D Ising model. A sustained increase of a SF with 
lattice size cannot be stronger than being proportional to 
V = L 3 , because an upper bound on each SF is obtained 
by setting all values in the sum of Eq. © equal to one. 

A. Equilibrium Results 

In this section we compile SF estimates from equilib- 
rium simulations of the 3D Ising and 3-state Potts model 
on L 3 lattices. Our simulations are carried out in a mul- 
ticanonical ensemble ^4|, covering a temperature range 
from /3 min = (infinite temperature) to /3 max > below 
the phase transition temperature of the respective sys- 
tem. Instead of relying on a recursion (see, e.g., ^2), the 



1. 3D Ising Model 

At the critical point the two-point function on an infi- 
nite lattice falls off with a power law, which defines the 
critical exponent ry: 

f{x) = < S (0) s(x) ) ~ \x r d+2 -" , \x\ - co . (14) 

This determines the low-momentum behavior of the 
Fourier transformation F(k). Namely, 

= J ^ e il2 \d-2+r, |-,|-«H-2-„ 

~ A"- 2 F , (fc) 

holds and, therefore, 

F{k)~\k\- b , b = 2-n, for |jfe|->0. (15) 

For fixed n we have k = 2i:n/ L and we find for any fixed 
value of n the finite size scaling (FSS) divergence 

F Pl - L b , 6 = 2-77, forL->co (16) 

with lattice size. 
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FIG. 1: SFs Fi, i — 1, . . . , 9 from simulations of the Ising FIG. 2: Finite size behavior of SF Fi from Ising model simu- 



model on an 80 lattice. 



lations on L lattices. 



The infinite volume phase transition temperature of 
the 3D Ising model is estimated to be (3 C — 0.22157(3), 
e.g., see ha. In our multicanonical simulations we cover 
the range from /3 m ; n = to /3 max = 0.25, well including 
the transition region. Table [Q gives an overview of the 
lattice sizes and the accumulated statistics (a sweep up- 
dates each spin once) together with our estimate of the 
maximum values F 1 max of the first SF, evaluated at the 
value (3 m . Error bars are given in parenthesis and apply 
to the last digits of the number in front. They are cal- 
culated with respect to a number of jackknife bins given 
by the first number in column two of the table, and the 
multicanonical re-weighting procedure uses the logarith- 
mic coding described in |17j . Three independent runs 
were carried out for the L = 80 lattice. Before starting 
with measurements we normally performed the number 
of sweeps of one measurement bin for reaching equilib- 
rium. This is sufficient because equilibration problems 
are mild in multicanonical simulations. Running time for 
each of our L = 80 simulations was about three months 
on a 2 GHz Athlon PC. The last column of table[Q gives 
the number of cycles 

(fie < Prnin) -> (Pe > Anax) -> (fie < Anin) , 

which the Markov process performed during the produc- 
tion run, where j3 e is the effective energy-dependent (3 of 
the multicanonical procedure. 

For our largest lattice the SFs 1-9 are plotted in Fig.Q^ 
where we restrict (3 to a neighborhood of the critical tem- 
perature. Each SF develops a clear peak, only that the 
peaks for the higher SFs are less pronounced than those 
for the lower. In particular the scale of the figure does not 
resolve the peaks for the SFs > 7 anymore. These peaks 
are found on a reduced scale and for each SF the FSS 
behavior l|16|) holds. However, the numerical accuracy 

decreases with increasing k . So we are content with 

simply analyzing the FSS behavior of SF 1. Fig.[2shows 
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b= 1.959 (12) 
Q = 0.82 
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FIG. 3: Fit of the F\ max maxima of table |T| to the FSS 
form (1161) . 



SF 1 for all our lattice sizes and the maxima values are 
collected in table Q] A two parameter fit to the form ll(i|) 
is shown in Fig.|3] It gives b = 1.959 (12) with a goodness 
of fit Q = 0.82 (for the definition of Q see, e.g., Ref. ^2), 
a result well compatible with the high precision estimates 
r\ = 0.0364 (5) given in the review article pj§| on critical 
phenomena and renormalization group theory. 



2. 3D 3-state Potts Model 

For the 3D 3-state Potts model one deals with a 
relatively weak first order phase transition at (3 C = 
0.2752720 (49), a value which averages two somewhat in- 
consistent (Q = 0.003 for the Gaussian difference test) es- 
timates of the literature (because of the inconsistency 
the error bars are averaged here and not reduced). In 
our multicanonical simulations we cover the range from 
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TABLE II: Statistics and SF maxima f\ max at f3 m for our 
equilibrium simulations of the 3D 3-state Potts model on L A 
lattices. 



L 


sweeps 


|jimax 




cycles 


20 


32 x 1.2 ■ 10 5 


19.00 (22) 


0.274273 


59 


30 


32 x 5.2 ■ 10 5 


38.11 (41) 


0.274924 


71 


40 


32 x 1.5 ■ 10 6 


60.30 (50) 


0.275116 


73 


50 


126 x 1.5 ■ 10 6 


80.46 (55) 


0.275181 


131 




0.274 0.2745 0.275 0.2755 0.276 0.2765 
P 

FIG. 4: SFs Fi, i — 1, . . . , 9 from simulations of the 3-state 
Potts model on a 50 3 lattice. 

Anin = to /?max = 0.33. An overview of the statistics 
and some results are given in table [n] similarly as before 
for the 3D Ising model in table [I] 

For our 50 3 lattice the SFs 1-9 are plotted in Fig. ^ 
where we restrict (3 to a neighborhood of the transition 
temperature. As for the 3D Ising model each SF develops 
a clear peak, but the shapes are significantly different. A 
relatively smooth increase is followed by a rather abrupt 
decrease. The lattice size dependence of SF 1 is depicted 
in Fig. [SJ which indicates (as expected) that the abrupt 
decrease develops into a discontinuity for L — > oo. The 
increase of the structure function maxima is irregular and 
smaller from L = 40 to L = 50 than from L = 30 to 
L = 40. Asymptotically for L — > oo a finite maximum 
value is expected in case of a first order phase transition. 
Within our limited lattice sizes this is not yet seen, but a 
power law fit (|16|) of the type of Fig. |3 which is the large 
L behavior in case of a second order transition, becomes 
entirely inconsistent: Q — 2.7 • 10~ n is the goodness of 
fit obtained. 



B. Quenches 

After outlining the equilibrium scenario, let us discuss 
the time evolution after a quench from the disordered 
into the ordered phase of the 3D 3-state Potts model. 




0.274 0.2745 0.275 0.2755 0.276 0.2765 
P 

FIG. 5: Finite size behavior of SF F\ from 3-state Potts model 
simulations on L 3 lattices. 

TABLE III: Repetitions of quenches from /3 = 0.2 to (3f for 
the 3D 3-state Potts model on I? lattices. 



P f \L: 


20 


40 


60 


80 


100 


120 


0.3 


640 


640 


640 


320 


320 


320 


0.27 


640 


320 


320 


32 







An overview of our statistics is given in table IIIII We 
quench from j3 = 0.2 to the /?/ value given in the table, 
which collects the numbers of repetitions of each quench. 
Error bars are then calculated with respect to 32 bins. 
Larger lattices exhibit self-averaging, so that one needs 
less repetitions than for smaller lattices. 

In previous work Q we have investigated the quench 
(3 = 0.2 — > 0.3 and its subsequent stochastic time evo- 
lution on lattices up to size 80 3 . Meanwhile we have 
extended the SF part of this investigation to lattices of 
size up to 120 3 and Fig. shows the time evolution of 
SF 1 after this quench. Note that we divide the SF by 
an extra volume factor in this figure. So its initial in- 
crease with lattice size is faster than ~ V , the maximum 
sustained increase we discussed in the first paragraph of 
section ITTTl For our largest lattices, L > 80, the increase 
appears to level off to precisely 

F x maJ£ (L) ~ V = L 3 for L^oo (17) 

Heuristically this behavior during spinodal decomposi- 
tion may be expected: The quench changes the tempera- 
ture in the entire lattice instantaneously. It is then plau- 
sible that the local contribution to the SF is, in the av- 
erage, everywhere the same. So one expects an increase 
~ V of the maxima. The initial overshooting may be 
explained by an increase of correlations with lattice size, 
which levels off once the lattice size exceeds the correla- 
tion length. 

To test how this growth of the signal proportional to 
the volume depends on the depth of the quench into the 
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FIG. 6: Time evolution of SF Fi for the 3-state Potts model 
on L 3 lattices after a quench from (3 — 0.2 — > 0.3. 



FIG. 8: Time positions of the SF F± maxima versus lattice 
size. 
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FIG. 7: Time evolution of SF Fi for the 3-state Potts model 
on L 3 lattices after a quench from (3 — 0.2 — > 0.28. 



ordered region, we performed a quench to a temperature 
closer to the transition temperature, (3 = 0.2 — > 0.28. As 
shown in Fig. E|we find the same phenomenon as before: 
The maximum sustained increase ~ V is initially over- 
shot. The growth of the signal is weaker than before, as is 
expected since the system does not change so drastically. 

In both figures we see that the time positions i max 
of the SF 1 maxima move towards larger values with 
increasing lattice size. For our two quenches i max (L) is 
plotted in Fig. |H1 on a log-log scale. With parameters 
ao and a\ both curves can consistently be fitted to the 
expected form 



1 {L) = a + a 1 L 2 



(18) 



As t is measured in units of sweeps, the number of spin 
updates per time unit does not depend on L. In spin sys- 
tems t is thus proportional to the physical time. After 
the quench into the ordered phase the infinite spin system 



cannot be equilibrated in any finite time, a fact known in 
condensed matter physics j2£|. The explanation for this 
phenomenon is that the systems grows initially compet- 
ing domains of three distinct orientations. To dissolve 
these domains by local random fluctuations until one of 
them dominates the entire lattice is a slow process, which 
requires of order L 2 time. 

Visualization of these domains faces difficulties, be- 
cause naive geometrical definitions do not work. Com- 
pare Fig. 8 of Ref. Q. For analogue Potts models this 
is overcome by the Fortuin-Kasteleyn [2l| cluster defini- 
tion, but there is no immediate generalization to gauge 
theories, although promising ideas have been published 
|22j. Here we do not investigate this question any fur- 
ther. We think that it is safe to assume that competing 
domains are in both, spin and gauge systems, the under- 
lying cause for the explosive growth of structure factors 
Fi (t) , which we encounter in their time evolution after a 
heating quench. 

Finally in this section, based on 640 repetitions Fig. [5] 
demonstrates that for a non-critical quench nothing more 
than a smooth transition from one equilibrium value to 
the next happens. Therefore the explosive growth of SFs 
is a unambiguous signal that /3/ is indeed in the ordered 
phase. 



IV. SU(3) 

We report results from quenches of pure SU(3) lattice 
gauge theory on N T N% lattices. Our statistics is sum- 
marized in tables IIVI and All quenches are from the 
initial value 6/g 2 = 5.5. The 4 x simulations of ta- 
ble \W\ were already reported in ^(j- The simulations for 
the other lattices are new. The difference between the 
tables is that for the lattices of table IIVI we follow the 
quench all the way to its equilibrium value at Tf, while 
for the lattices of table [V] we calculated only the initial 
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FIG. 9: Time evolution of SFs F h i = 1, ... ,4 for the 3-state 
Potts model on a 40 3 lattice after a (non-critical) quench from 
(5 = 0.2 -> 0.27. 



TABLE IV: Quenches from 6/# 2 = 5.5 to 6/g 2 for pure SU(3) 
lattice gauge theory (n denotes the number of repetitions). 



Lattice 


Tt/T e 


6/0? 


n 


T f /T c 


6/3? 


n 


4 x 16 3 


1.250 


5.802740 


10 000 


1.568 


5.920000 


10 000 


4 x 32 3 








1.568 


5.920000 


4 000 


4 x 64 3 








1.568 


5.920000 


170 


6 x 24 3 


1.250 


6.022334 


6 000 


1.568 


6.165427 


6 000 


8 x 32 3 


1.250 


6.206036 


3 000 


1.568 


6.364572 


3 000 



TABLE V: Initial quenches from 6/g 2 = 5.5 to 6/<?? for pure 
SU(3) lattice gauge theory (n as in table ITVl . 



Lattice 




6/<7? 


n 


T f /T c 


6/P? 


n 


4 x 32 3 


1.250 


5.802740 


3 000 








4 x 64 3 


1.250 


5.802740 


140 








6 x 48 3 


1.250 


6.022334 


600 


1.568 


6.165427 


750 


6 x 60 3 


1.250 


6.022334 


200 
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6.165427 


200 


8 x 56 3 


1.250 


6.206036 
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FIG. 10: Time evolution of SF Fi/V for SU(3) lattice gauge 
theory on 4 x iV 3 lattices after a quench 6/g 2 = 5.5 — > 5.92. 



increase of the SFs as needed for the determinations of 
critical modes in section HV Bl 

The new data serve to study the quantum continuum 
limit a — > (in physical units like fermi) . The final values 
g 2 - of the bare coupling constants are chosen, so that the 
values of Tf/T c stay at the fixed ratios given in the table. 
For this we take (substantial) corrections to the two-loop 
equation of Lambda lattice into account, which follow 
from renormalization group results tabulated in Ref. |23| . 
As the use of tables is tedious, we like to mention that 
with an accuracy of 0.5% and better our Tf/T c values 
are reproduced by using the formula 



A L ( 5 2 ) = A£V)A( 5 2 ) 
where Af^g 2 ) is given by (e.g., pij) 



(19) 



For N T = 4, fixed, Fig. shows the divergence of 
the SF 1 maxima with increasing lattice size N% as well 
as a {"""'(Na) ~ N 2 behavior in complete analogy to 
our results for Potts models. All the lattices of Fig. EH 
are quenched to the bare coupling constant gj = 6/5.92. 
Therefore the time scale of the Markov process (deter- 
mined by the Boltzmann factors) is the same on all these 
lattices and up to an unknown multiplicative factor iden- 
tified with that of a dissipative, non-relativistic dynam- 
ics. Non-relativistic does not necessarily mean that the 
propagation of the signal through the lattice is slow. In 
the contrary, Galilee transformations set no upper limit 
on speeds. Our quench changes the temperature instan- 
taneously through the entire lattice, while the subsequent 
propagation of the response proceeds through local inter- 
actions. 



Ar 



2^-&i/(MS) e -l/(2fc 0S 2 ) 



Finite Volume Continuum Limit 



with bo = ^- < - 



h — 34 

°1 - T V 16-7T 2 



and 



A( 5 2 ) = l + a ie ~ a */° 2 +a 3 g 2 + a 4 g 4 



with ai = 71553750, 
and a 4 = 0.5089052. 



a 2 



19.48099 a 3 = -0.03772473, 



In the following we illustrate the approach of the limit 
a — > 0, L = constant, Tf/T c = constant, by increasing 
N T from 4 to 6 to 8 and the volume N% from N a = 16 
to 24 to 32, so that the ratio N a /N T stays constant. 
The ratio of temperatures Tf/T c is kept constant by us- 
ing Eq. i|19fl to determine the appropriate bare coupling 
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FIG. 11: Time evolution of SF Fi/Fi,/ for SU(3) lattice gauge 
theory on N r N% lattices of constant physical volume of a 
quench to Tf /T c = 1.25. 



10 
8 
6 
4 
2 




4x1 6„ 




200 400 

t' 



600 



FIG. 12: Time evolution of SF Fi/Fij for SU(3) lattice gauge 
theory on N T N% lattices of constant physical volume of a 
quench to T f /T c = 1.568. 



constants values for each N T . Due to the divergence of 
(bare) Polyakov loop correlations we face a renormaliza- 
tion problem, which we overcome by dividing all SFs i*j 
by their equilibrium values at Tf, Fij. The time-scale 
situation changes too, because we have to use different 
bare coupling constants values for different N T . As one 
knows that a finite physical volume equilibrates in a fi- 
nite time, we fix this normalization problem by rescaling 
the time axis to 



t' 



t 



\ t (N T ,T f /T c ) 



(20) 



so that all maxima fall on top of one another. We do not 
lose information as we anyhow do not know the overall 
normalization factor for our time scale. 

Figures El an d 113 show the time evolution of the 
Fi/Fij SFs for our two Tf/T c values. The time axis 
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FIG. 13: SU(3) determination of p c for T f /T c = 1.25. 



of our original measurements in units of sweeps are re- 
lated to those used in Fig.^Jby the Xt(N T , 1.25) factors 
1 : 2.655 : 5.457 for N T the values 4:6:8, respec- 
tively. For Fig.^Jthe corresponding Xt(N T , 1.568) ratios 
are 1 : 2.768 : 6.362. The maxima of the curves decrease 
when increasing N T from 4 to 6 to 8. As the decrease 
slows down with increasing lattice size, there is some ev- 
idence for an approach to a shape, which represents the 
continuum limit. 



B. Debye Screening Mass 

The current understanding of the early time evolution 
of systems out of equilibrium is largely based on inves- 
tigating stochastic equations which are dynamical (time 
dependent) generalizations of the Landau-Ginzburg ef- 
fective action models of the static (equilibrium) the- 
ory UIS]. For model A the linear approximation results 
in the following equation for a SF: 



dF(p,t) 



dt 



2w(p)F(p,t) 



with the solution 



F(p,t) = F(p,t = 0)exp(2u;(p)t) 
w(p) > for |p| > p c 



(21) 



(22) 



where p c > is a critical momentum. Originally the 
linear theory was developed for model B |26l |27| . Details 
for model A can be found in Ref . 0, • 

From our measurements of F(p, t) on the N T = 4, 
6 and 8 lattices we find straight line fits to the form 
oj(p) = clq +aip 2 , p — \p\ with a negative slope a%. They 
determine the critical momentum p c as the value where 
oj{p) changes its sign. The fits for Tf/T c = 1.25 are 
shown in Fig.USand for T f /T c = 1.568 in Fig-ITU where 
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FIG. 14: SU(3) determination of p c for T f /T c = 1.568. 



TABLE VI: Fit results for p c /T c . 



Lattice size: 


N T = 4 


N T = 6 


AV = 8 


oo 


T f /T c = 1.25 : 


1.613 (18) 


1.424 (26) 


0.87 (26) 


1.023 (85) 


T f /T c = 1.568 : 


2.098 (19) 


2.095 (22) 


2.40 (12) 


2.038 (73) 



we introduced 

w'(p) = \ t (N r ,T f /T c )u(p) . (23) 

This definition absorbs the shift 120(1 of the time scale, so 
that u)'(p) t' — u)(p)t holds. It is only in the primed vari- 
ables that one realizes an eventual approach to the con- 
tinuum limit from Figs. ll3l and ll4l In particular note that 
for Tf/T c = 1.568 the N T = 6 and 8 fits are within statis- 
tical errors identical. The obtained values for p c (N r )/T c 
are listed in table IVII The (finite volume) continuum 
limit is extrapolated by fitting these values to the form 

Pc(N T ) pc_ const 

with the results given in the last column of table IVII 

Relying on a phenomenological analysis by Miller and 
Ogilvie 9], p c is related by 

Vs Pc (25) 

to the Debye screening mass at the final temperature Tf 
after the quench. We get 

m D = 1.77 (15) T c for T f /T c = 1.25, (26) 
m D = 3.53 (13) T c for T f /T c = 1.568. (27) 

The value at Tf/T c — 1.568 is in excellent agreement with 
a determination of mu (T) from a best-fit analysis of the 
large distance part of the color singlet free energies |2flj . 
This supports that the simulated dynamics bears physi- 
cal content. Our estimate at Tj/T c = 1.25 is by a factor 



of two smaller than the one of Ref. [2^ . This is not really 
a surprise, because Tj JT C = 1.25 is close to the spinodal 
cndpoint, so that the derivation [j| of the relationship 
((25(1 is no longer valid. 

For pure SU(3) lattice gauge theory T c = 265 (1) MeV 
holds, assuming a — 420 MeV for the string tension, 
while for QCD the cross-over temperature appears to be 
around T c « 165 MeV, see Ref. [3(j for reviews. Us- 
ing for simplicity T c = 200 MeV to illustrate the mag- 
nitudes, the temporal lattice size is then about 1 fermi 
at T c . The spatial sizes of our lattices used in this sec- 
tion reach up to (8 fermi) 3 . At the Tj values the edge 
lengths are shortened by the corresponding T c /Tt fac- 
tors. I.e., the volume is (6.4fermi) 3 for Tf/T c = 1.25 and 
(5.10 fermi) 3 for T f /T c = 1.568. The screening length 
associated with the Debye mass, £d = I /mo, is then 
approximately 0.6 fermi at Tf/T c = 1.25 and 0.3 fermi 
at T f /T c = 1.568. The illust ration of the finite vol- 
ume continuum limit in section IIV Al was for lattices of 
size (4fermi) 3 at T c , i.e., (3.2 fermi) 3 at T f /T c = 1.25 
and (2.55 fermi) 3 at Tt/T c — 1.568. Our volumes are 
smaller than the envisioned deconfined region of about 
(10 fermi) 3 in relativistic heavy ion experiments. Due to 
periodic boundary conditions one may expect that MC 
simulations on smaller lattices are representative for the 
central region of the larger volume. Our result is that 
the Debye screening length is short on the scale of the 
deconfined region. 

C. Measurements near Structure Factor Maxima 
versus Deconfined Equilibrium 

For SU(3) gauge theory the triality of Polyakov loops 
with respect to the Z 3 center of the gauge group takes 
the place of three distinct spin orientation. Although 
a satisfactory cluster definition does not exist for gauge 
theories, the underlying mechanism of competing vacuum 
domains is expected to be similar as in the spin models. 
To study their influence on Polyakov loop correlations 
and on the gluonic energy e and pressure p densities, we 
calculate these quantities at times t < i max as well as at 

Our structure function measurements gave "on the fly" 
two-point correlations between Polyakov loops defined by 

C (d, t) = (P(0, t)P(d, t)) L - «|P(0, t)\) L f (28) 

where the averaging procedures are those we discussed 
after Eq. (0 and d = 1,2, y/2, .... The value of these 
results is somewhat limited, because our focus was not 
on good equilibrium results and the stored data do not 
allow to project onto particular channels of the free en- 
ergy of static quarks (which lead to larger correlation 
lengths than those obtained) . For several values of d we 
plot in Fig. ^] the time development of C (d) on our 
largest lattice using the t' time scale ((20(1 . The corre- 
lations assume maxima at about the same time values 
imax for which the SFs peak, although less pronounced. 
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FIG. 15: Time dependence of Polyakov loop correlations 
C (d) on the 8 x 32 3 lattice for T f /T c = 1.25 and various 
d values. 
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FIG. 16: Fall-off behavior of the Polyakov loop correlations 
of Fig. 1151 at different times. 



In Fig. ^j] we plot the (i-dependence for the time val- 
ues 0.5t max , i max and 5i max . At 5t max fits of the form 
C (d) ~ exp(— mpad)/(ad), where a is the lattice spac- 
ing, give the mp estimates which are collected in ta- 
ble IVIII The last column of this table gives infinite vol- 
ume estimates obtained from fits of the form i|24|) . In 
contrast to that large correlations are found at 0.5 i max 
and t max , which are fully consistent with a power law. 

The equilibrium procedure for calculating the glu- 
onic energy e and pressure p densities is summarized in 
Ref. [23l l3l| (in earlier work the pressure exhibited 
a non-physical behavior after the deconfining transition 
and the energy density approached the ideal gas limit 
too quickly because the anisotropy coefficients were cal- 
culated perturbatively). We denote expectation values of 
space-like plaquettes by P a and those involving one time 



TABLE VII: Fit results for m P /T c at 5t n 



Lattice size: 


4 x 16 3 


6 x 24 3 


8 x 32 3 


00 


Tf/T c = 1.25 : 


3.27 (19) 


3.70 (20) 


4.43 (27) 


5.23 (49) 


T f /T c = 1.568 : 


4.82 (61) 


5.33 (35) 


6346 (38) 


7.70 (95) 
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FIG. 17: SU(3) gluonic energy density on the 4 x 16 3 lattice: 
(a) with Po calculated from the time series after the quench 
and (b) using equilibrium values for Po. 



link by P T . The energy density and pressure can then be 
cast into the form 



r 



e + p _ 8N C N? 

and 
e — 'ip 



1 ~ yM a ) ~ c T (a)] 



(P CT -P T ) (29) 



y4 



= 12iV c iV 4 Ma) - c T (a)} [2P - (Pa + P T )} , 



(30) 

where Po is the plaquette expectation value on a sym- 
metric (T = 0) lattice, and the anisotropy coefficients 
c <7iT (a) are defined by: 



Ca,r(a) 



2 

<J,T 



(31) 



They are related to the QCD /3-function and can be 
calculated using Pade fits of j2j|. To normalize to zero 
temperature, plaquette values from the symmetric N T = 
N a lattice are needed in Eq. (|3Tj|) . As one stays within the 
confined phase on the symmetric lattice its equilibration 
after the quench is fast. Therefore it is enough to use 
equilibrium values of Po at fif after the quench. This is 
illustrated in Fig. H3 

In Fig. 1181 we show the time evolution of the gluonic 
energy densities (upper curves) and pressure densities 
(lower curves) for the Tf/T c — 1.25 quench on our 4x 16 3 , 
6 x 24 3 and 8 x 32 3 lattices using the rescaled time defini- 
tion H20|) . The curves for the last two lattices fall almost 



11 











e/T 4 , N T =4 
e/T 4 , N T =6 










e/T 4 , N T =8 
3p/T 4 , N T =4 










3p/T 4 , N T =8 - 



200 300 400 500 600 

f 

FIG. 18: SU(3) gluonic energy densities and pressures at 
T f /T c = 1.25. 

on top of one another, indicating their neighborhood to 
the continuum limit. The approach to the final equi- 
librium values is rather smooth. Gluonic energy density 
mean values at i max are less than 1/4 of their final values, 
while the pressure density is at about 1/3. In contrast 
to the shift in the mean value, the width of the distribu- 
tions are almost the same at i max and in the deconfincd 
equilibrium. Results for the Tf/T c — 1.568 quench are 
quite similar. 



V. SUMMARY AND CONCLUSIONS 

In equilibrium at temperatures much higher than the 
deconfinement temperature T c the perturbative prescrip- 
tion of QCD is that of a weakly coupled gas of quasi- 
particles. In contrast to that recent experiments at the 
BNL relativistic heavy-ion collider (RHIC) show coher- 
ence in particle production and strong collective phenom- 
ena, which are well described by the model of a near- 
perfect, strongly coupled fluid |33j |. Non-perturbative ef- 
fects are expected to play some role in the prescription of 
equilibrium QCD at temperatures reached at the RHIC. 
For the T//T c = 1.25 and T f /T c = 1.568 temperatures 
investigated in this paper equilibrium lattice calculations 



indicate indeed corrections (compare Fig. 7 of [23]). How- 
ever correlations are typically over ranges much smaller 
than the deconfined region, compare our estimates of the 
Debye screening mass m£>(T/). The agreemet of our md 
value at Tf /T c — 1.568 with direct equilibrium estimates 
[2^ | give confidence that model A dynamics reflects phys- 
ical features. 

If the phenomenological description of a strongly cou- 
pled plasma implies correlations over distances exceed- 
ing one fermi, the time evolution of our structure factors 
(SFs) depicted in Fig. IrJl 171 ISI and ITU1 suggest a scenario in 
which the deconfined equilibrium phase has actually not 
been reached at the RHIC, but the heating process gets 
stuck during the time period of explosive growth of the 
SFs. While this explains correlation over distances much 
larger than one fermi, it also provides an unambiguous 
signal for the existence of the deconfining phase: Fig. [5] 
demonstrates that the explosive growth is absent for a 
non-critical quench. 

In real QCD there are two effects which prevent the 
divergence of the equilibration time shown in Fig. [5] 
(1) Quarks break the Z 3 symmetry of the SU(3) gauge 
group, similarly as a magnetic field breaks the degener- 
acy of the spins in the 3D 3-state Potts model. The final 
magnitude of the equilibration time depends then on the 
strength of the breaking as illustrated in Ref. 2] for a 
weak magnetic field. (2) At the RHIC the physical vol- 
ume is finite, so that even in case of an exact symmetry 
the equilibration time is finite. So the scenario that the 
system gets stuck during the spinodal decomposition of 
its vacuum structure could only be based on phenomeno- 
logical observations. Questions like how a perfect fluid 
may look during the period of spinodal decomposition 
arise. Minkowski space simulations of hyperbolic dif- 
ferential equations, which emerge from effective actions 
for Polyakov loops 0, |35| , may shed light on the ques- 
tion whether features observed in the paper are special 
to Glauber dynamics or of some universal nature. 
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